home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 20 / Cream of the Crop 20 (Terry Blount) (1996).iso / math / alged34.zip / ALGEDSRC.ZIP / ALGPOLY.C < prev    next >
C/C++ Source or Header  |  1996-06-06  |  16KB  |  540 lines
  1. /*--------------------------------------------------------------------
  2.    Alged:  Algebra Editor    henckel@vnet.ibm.com
  3.  
  4.    Copyright (c) 1994 John Henckel
  5.    Permission to use, copy, modify, distribute and sell this software
  6.    and its documentation for any purpose is hereby granted without fee,
  7.    provided that the above copyright notice appear in all copies.
  8. */
  9. #include "alged.h"
  10. /*-----------------------------------------------------------------
  11.    make-tree convert a table of coefficients to a node tree.
  12.    Note
  13.    the coefficients are REUSED (NOT copied), so don't free them.
  14. */
  15. node *maketree(node **coef, int sz, node *base) {
  16.   int i;
  17.   node *p,*tmp;
  18.  
  19.   p = coef[0];
  20.   for (i=1; i<sz; ++i) {
  21.     if (coef[i]->kind==NUM && !coef[i]->value) {
  22.       freenode(coef[i]);             /* throw away zeros */
  23.     }
  24.     else {
  25.       tmp = p;
  26.       p = newoper(ADD);
  27.       p->lf = tmp;
  28.       p->rt = tmp = newoper(MUL);
  29.       tmp->lf = coef[i];             /* add coef[i]*base^i */
  30.       if (i>1) {
  31.         tmp->rt = newoper(EXP);
  32.         tmp->rt->lf = deepcopy(base);
  33.         tmp->rt->rt = newnum(i);
  34.       }
  35.       else tmp->rt = deepcopy(base);
  36.     }
  37.   }
  38.   return p;
  39. }
  40.  
  41. /*-----------------------------------------------------------------
  42.    looks inside 'a' for any expression b or b^N.
  43.    if found, it changes it to a 1 (ONE) and returns the exponent.
  44.    Note, if N is not an integer (it's a fraction or expression)
  45.    then findbase ignores it.
  46.    When findbase returns 0, then  you can assume 'a' is unchanged.
  47. */
  48. int findbase(node *a, node *b) {
  49.   int i,x=1;
  50.  
  51.   if (equal(a,b) ||                         /* any expression */
  52.       a->kind==EXP && a->rt->kind==NUM &&       /* EXPONENTS */
  53.       (x = a->rt->value, whole(x)) &&
  54.       x >=0 && equal(a->lf,b)) {
  55.     movenode(a,newnum(1));
  56.     return x;
  57.   }
  58.   if (a->kind==MUL) {                     /* recurse on MULTIPLY */
  59.     x = findbase(a->lf,b);
  60.     if (x) return x;
  61.     x = findbase(a->rt,b);
  62.     return x;
  63.   }
  64.   return 0;
  65. }
  66.  
  67. /*-----------------------------------------------------------------
  68.    add entry to coef table          grow and init table
  69. */
  70. #define addcoef(k,a) do {                 \
  71.   node *tmp;                              \
  72.   if (sz<=i) {                            \
  73.     coef = realloc((void*)coef,(i+1)*sizeof*coef);   \
  74.     for ( ; sz<=i; ++sz)                  \
  75.       coef[sz] = newnum(0);               \
  76.   }                                       \
  77.   tmp = coef[i];                          \
  78.   coef[i] = newoper(k);                   \
  79.   coef[i]->lf = tmp;                      \
  80.   coef[i]->rt = a;                        \
  81. } while(0)
  82.  
  83. /*-----------------------------------------------------------------
  84.    MAKETABLE convert a tree to a table of coefficients which is indexed
  85.    by the power of a given base.  For example, a*x^2 + 3*x - b would be
  86.       0 -> 0 - b
  87.       1 -> 0 + 3
  88.       2 -> 0 + a
  89.    the size returned is the index of the last non-zero coefficient + 1
  90.    Note:
  91.    the tree, p, is expected to have correct association.
  92.    the coefficients always have a leading zero.
  93.    the tree is not altered or referred to by the table.
  94.    on return, you are guaranteed that result[i]->kind is one of NUM,
  95.       ADD, SUB and if it is NUM then it is zero.
  96. */
  97. node **maketable(node *base, node *p, int *size) {
  98.   node *a;
  99.   int i,sz;
  100.   node **coef;
  101.  
  102.   coef = malloc(sizeof*coef);        /* initial size = 1 */
  103.   checknull(coef);
  104.   *coef = newnum(0);
  105.   sz = 1;
  106.   a = deepcopy(p);                /* make a working copy */
  107.   while (1) {
  108.     i = findbase(a,base);
  109.     if (!i && aop(a->kind)) {
  110.       i = findbase(a->rt,base);
  111.       addcoef(a->kind,a->rt);
  112.     }
  113.     else {
  114.       addcoef(ADD,a);
  115.       break;
  116.     }
  117.     a = a->lf;
  118.   }
  119.   *size = sz;
  120.   return coef;
  121. }
  122.  
  123. /*-----------------------------------------------------------------
  124.    polynomial long division
  125.  
  126.    base is the expression on which the division is done (e.g. x)
  127.    nm and dn are numerator and denominator
  128.    returns:
  129.    a new tree if success, else returns NULL
  130. */
  131. node *polydiv(node *base, node *nm, node *dn) {
  132.   node **num,**den,**quo,**rem;
  133.   int szn,szd,szq,szr;
  134.   int i,j;
  135.   node *tmp,*p;
  136.  
  137.   num = maketable(base,nm,&szn);
  138.   den = maketable(base,dn,&szd);
  139.  
  140.   /*  If degree of numerator is less than denominator, or
  141.       if the denominator does not contain base, then return failure */
  142.  
  143.   if (szn < szd || szd < 2) {
  144.     for (i=0; i<szn; ++i) freetree(num[i]);
  145.     for (i=0; i<szd; ++i) freetree(den[i]);
  146.     free(num); free(den);
  147.     return NULL;
  148.   }
  149.   szq = szn-szd+1;
  150.   quo = malloc(szq*sizeof*quo);
  151.   checknull(quo);
  152.   szr = szd-1;
  153.   rem = malloc(szr*sizeof*rem);
  154.   checknull(rem);
  155.  
  156.   /*  Divide every coefficient by leading coeff in denominator */
  157.  
  158.   p = den[szd-1];
  159.   if (!(p->kind==ADD && p->rt->kind==NUM && p->rt->value==1.0)) {
  160.     for (i=0; i<szn; ++i) {
  161.       if (num[i]->kind != NUM) {      twirl();
  162.         tmp = num[i];
  163.         num[i] = newoper(DIV);
  164.         num[i]->lf = tmp;
  165.         num[i]->rt = deepcopy(p);
  166.       }
  167.     }
  168.     for (i=0; i<szd-1; ++i) {
  169.       if (den[i]->kind != NUM) {      twirl();
  170.         tmp = den[i];
  171.         den[i] = newoper(DIV);
  172.         den[i]->lf = tmp;
  173.         den[i]->rt = deepcopy(p);
  174.       }
  175.     }
  176.   }
  177.  
  178.   /*  Construct the coefficients of the quotient */
  179.  
  180.   for (i=1; i<=szq; ++i) {
  181.     quo[szq-i] = deepcopy(num[szn-i]);
  182.     for (j=1; j<i; ++j) if (i-j < szd) {      twirl();
  183.       tmp = quo[szq-i];
  184.       quo[szq-i] = newoper(SUB);
  185.       quo[szq-i]->lf = tmp;
  186.       quo[szq-i]->rt = tmp = newoper(MUL);
  187.       tmp->lf = deepcopy(den[szd-i+j-1]);
  188.       tmp->rt = deepcopy(quo[szq-j]);
  189.     }
  190.   }
  191.  
  192.   /*  Construct the coefficients of the remainder */
  193.  
  194.   for (i=0; i<szr; ++i) {
  195.     rem[i] = deepcopy(num[i]);
  196.     for (j=0; j<=i; ++j) if (i-j < szq) {      twirl();
  197.       tmp = rem[i];
  198.       rem[i] = newoper(SUB);
  199.       rem[i]->lf = tmp;
  200.       rem[i]->rt = tmp = newoper(MUL);
  201.       tmp->lf = deepcopy(den[j]);
  202.       tmp->rt = deepcopy(quo[i-j]);
  203.     }
  204.   }
  205.  
  206.   /*  Convert the quo and rem tables into a node tree */
  207.  
  208.   p = newoper(ADD);
  209.   p->rt = newoper(DIV);
  210.   p->rt->lf = maketree(rem,szr,base);
  211.   p->rt->rt = deepcopy(dn);
  212.   p->lf = maketree(quo,szq,base);
  213.  
  214.   /*  Calculate remainder trivial stuff */
  215.  
  216.   while (calcnode(p->rt,1));
  217.  
  218.   /*  Clean up and return */
  219.  
  220.   for (i=0; i<szn; ++i) freetree(num[i]);
  221.   for (i=0; i<szd; ++i) freetree(den[i]);
  222.   free(num);
  223.   free(quo);
  224.   free(rem);
  225.   free(den);
  226.   return p;
  227. }
  228.  
  229. /*-----------------------------------------------------------------
  230.    polycoef - collect the coefficients of a polynomial
  231. */
  232. void polycoef(node *base, node *p) {
  233.   node **coef;
  234.   int sz;
  235.   int i;
  236.  
  237.   if (!p || !base || p->kind==EQU || base->kind==EQU) return;
  238.  
  239.   coef = maketable(base,p,&sz);
  240.   if (sz < 2) {
  241.     if (sz) freetree(coef[0]);
  242.     free(coef);
  243.     return;
  244.   }
  245.   for (i=0; i<sz; ++i)              /* calculate the coefficients */
  246.     while (calcnode(coef[i],1));
  247.  
  248.   movenode(p,maketree(coef,sz,base));     /* replace p */
  249.   free(coef);
  250. }
  251.  
  252. /*-----------------------------------------------------------------
  253.    use quadratic equation to factor a degree-2 polynomial
  254.                                      2                       2
  255.      2                     b + sqrt(b - 4ac)       b - sqrt(b - 4ac)
  256.    ax + bx + c  ==>  a(x + -----------------) (x + -----------------)
  257.                                  2a                       2a
  258. */
  259. void quadratic(node *base,node *p) {
  260.   node **coef,*a;
  261.   int sz;
  262.   int i;
  263.  
  264.   coef = maketable(base,p,&sz);
  265.  
  266.   if (sz != 3) {
  267.     for (i=0; i<sz; ++i) freetree(coef[i]);
  268.     free(coef);
  269.     return;
  270.   }
  271.  
  272.   /* Construct the two binomials p = (ax + u)(x + v)  */
  273.  
  274.   a = newoper(ADD);
  275.   a->lf = deepcopy(coef[1]);                     /* b */
  276.   a->rt = newoper(EXP);
  277.   a->rt->lf = newoper(SUB);
  278.   a->rt->lf->lf = newoper(EXP);                  /* b^2 */
  279.   a->rt->lf->lf->lf = deepcopy(coef[1]);
  280.   a->rt->lf->lf->rt = newnum(2);
  281.   a->rt->lf->rt = newoper(MUL);